Can we predict one’s biological age?

Introduction to Linear Modelling

Libraries

library(here)       # easy directory management, as here("YOUR/FILE/PATH")
library(data.table) # fast data wrangling
library(magrittr)   # use of pipe operator `%>%`
library(ggplot2)    # graphics
library(nortest)    # normality assessment test

Load data

We will be analysing a subset of the whole blood data from 656 human individuals, made available by Hannun et al. (2013).

# base R
mdna <- read.csv(here("data/methDNA_data_Hannun_etal_2013.csv"), header = TRUE)

# alternative: using data.table's fread()
# mdna <- fread(here("data/methDNA_data_Hannun_etal_2013.csv"), header = TRUE)

# look at first 8 cols from the data
head(mdna[, 1:8])

##      GEO.ID Age Gender Ethnicity cg17568996 cg02245378 cg14093936 cg15642326
## 1 GSM989827  67      F Caucasian 0.08144546 0.08951550  0.3297176  0.2666367
## 2 GSM989828  89      F Caucasian 0.06307138 0.09956124  0.2693837  0.2525504
## 3 GSM989829  66      F Caucasian 0.12086410 0.10301870  0.3216346  0.2430219
## 4 GSM989830  64      F Caucasian 0.09724835 0.08985170  0.3458824  0.2730094
## 5 GSM989831  62      F Caucasian 0.12239770 0.09564840  0.3497202  0.2836436
## 6 GSM989832  87      F Caucasian 0.19961670 0.16343000  0.3575867  0.2096256

Variables:

• GEO.ID: Unique identifier for each one of the 656 participants
• Age: Age of participants, in years
• Gender: Gender of each individual, binary values for female (F) and male (M)
• Ethnicity: Ethnic group of each individual, binary identifying European participants (Caucasian) and Mexican participants (Hispanic)
• cgXXXXXXXX: CpG marker, each one indicating the methylation values, in percentage

# get a vector with all column names with their pattern starting with "cg" (i.e., the names of all CpG probes)
vec_probe <- grep("^cg", colnames(mdna), value = TRUE)

Quick EDA

Things to look for:

1. Look at the population-related variables
   • Are there sampling differences between gender and/or ethnicity?
2. Look at age distribution
   • Does age change with gender/ethnicity?
3. Look at methylation levels
   • Should we transform the data?
   • Are the values correlated with age?
   • Could these values by affected by gender/ethnicity?

Demographics and age

# Look through the group categories
table(mdna$Gender)

## 
##   F   M 
## 338 318

table(mdna$Ethnicity)

## 
## Caucasian  Hispanic 
##       426       230

# gender and ethnicity
table(mdna$Ethnicity, mdna$Gender)

##            
##               F   M
##   Caucasian 222 204
##   Hispanic  116 114

# proportions
round(prop.table(table(mdna$Ethnicity, mdna$Gender), margin = 2), 2)

##            
##                F    M
##   Caucasian 0.66 0.64
##   Hispanic  0.34 0.36

round(prop.table(table(mdna$Ethnicity, mdna$Gender), margin = 1), 2)

##            
##                F    M
##   Caucasian 0.52 0.48
##   Hispanic  0.50 0.50

# Pearson's Chi-squared test
chisq.test(table(mdna$Gender, mdna$Ethnicity))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(mdna$Gender, mdna$Ethnicity)
## X-squared = 0.10788, df = 1, p-value = 0.7426

# Gender and Age
boxplot(mdna$Age ~ mdna$Gender)

quantile(subset(mdna, Gender == "F")$Age)

##   0%  25%  50%  75% 100% 
##   21   56   64   75  101

quantile(subset(mdna, Gender == "M")$Age)

##    0%   25%   50%   75%  100% 
## 19.00 52.00 65.00 75.75 96.00

# Mann-Whitney U test
wilcox.test(mdna$Age ~ mdna$Gender)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mdna$Age by mdna$Gender
## W = 55605, p-value = 0.4425
## alternative hypothesis: true location shift is not equal to 0

# ggplot2 alternatives

# boxplot:
ggplot(mdna, aes(x = Gender, y = Age)) +
  geom_boxplot()

# density plot:
ggplot(mdna, aes(Age, fill = Gender)) +
  geom_density(alpha = 0.5)

# Ethnicity and Age
boxplot(mdna$Age ~ mdna$Ethnicity)

quantile(subset(mdna, Ethnicity == "Caucasian")$Age)

##   0%  25%  50%  75% 100% 
##   19   61   71   79  101

quantile(subset(mdna, Ethnicity == "Hispanic")$Age)

##   0%  25%  50%  75% 100% 
##   23   47   56   62   84

# Mann-Whitney U test
wilcox.test(mdna$Age ~ mdna$Ethnicity)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mdna$Age by mdna$Ethnicity
## W = 79004, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

# ggplot2 alternatives

# boxplot:
ggplot(mdna, aes(x = Ethnicity, y = Age)) +
  geom_boxplot()

# density plot:
ggplot(mdna, aes(Age, fill = Ethnicity)) +
  geom_density(alpha = 0.5)

DNA methylation values

# base summary
summary(mdna[, vec_probe])

##    cg17568996        cg02245378        cg14093936       cg15642326    
##  Min.   :0.03878   Min.   :0.06202   Min.   :0.2360   Min.   :0.1110  
##  1st Qu.:0.09808   1st Qu.:0.09436   1st Qu.:0.3190   1st Qu.:0.2419  
##  Median :0.12126   Median :0.10358   Median :0.3444   Median :0.2601  
##  Mean   :0.12249   Mean   :0.10704   Mean   :0.3436   Mean   :0.2596  
##  3rd Qu.:0.14284   3rd Qu.:0.11700   3rd Qu.:0.3682   3rd Qu.:0.2795  
##  Max.   :0.27857   Max.   :0.26465   Max.   :0.5108   Max.   :0.4882  
##    cg24991452        cg04581377         cg01820374       cg06898923     
##  Min.   :0.01363   Min.   :0.002807   Min.   :0.1625   Min.   :0.02812  
##  1st Qu.:0.03247   1st Qu.:0.039772   1st Qu.:0.2537   1st Qu.:0.05137  
##  Median :0.04142   Median :0.044587   Median :0.2796   Median :0.05834  
##  Mean   :0.04235   Mean   :0.044973   Mean   :0.2794   Mean   :0.05933  
##  3rd Qu.:0.05175   3rd Qu.:0.050517   3rd Qu.:0.3049   3rd Qu.:0.06635  
##  Max.   :0.07500   Max.   :0.090808   Max.   :0.4095   Max.   :0.18484  
##    cg11579069        cg19560758       cg17385448       cg03271907    
##  Min.   :0.03447   Min.   :0.1396   Min.   :0.1857   Min.   :0.8770  
##  1st Qu.:0.06001   1st Qu.:0.2275   1st Qu.:0.3527   1st Qu.:0.9215  
##  Median :0.06579   Median :0.2463   Median :0.3789   Median :0.9335  
##  Mean   :0.06688   Mean   :0.2487   Mean   :0.3823   Mean   :0.9329  
##  3rd Qu.:0.07175   3rd Qu.:0.2685   3rd Qu.:0.4128   3rd Qu.:0.9447  
##  Max.   :0.10725   Max.   :0.4399   Max.   :0.5631   Max.   :0.9933  
##    cg04081402       cg26614073       cg17895149       cg22736354    
##  Min.   :0.7150   Min.   :0.2468   Min.   :0.1666   Min.   :0.1483  
##  1st Qu.:0.8334   1st Qu.:0.4023   1st Qu.:0.4865   1st Qu.:0.2878  
##  Median :0.8524   Median :0.4361   Median :0.5050   Median :0.3140  
##  Mean   :0.8504   Mean   :0.4337   Mean   :0.5043   Mean   :0.3130  
##  3rd Qu.:0.8705   3rd Qu.:0.4669   3rd Qu.:0.5247   3rd Qu.:0.3381  
##  Max.   :0.9591   Max.   :0.6279   Max.   :0.7508   Max.   :0.5102  
##    cg07549715       cg17592812        cg18463417       cg26832211    
##  Min.   :0.2752   Min.   :0.01744   Min.   :0.7925   Min.   :0.7379  
##  1st Qu.:0.4750   1st Qu.:0.04622   1st Qu.:0.8500   1st Qu.:0.8199  
##  Median :0.4998   Median :0.05143   Median :0.8593   Median :0.8430  
##  Mean   :0.5007   Mean   :0.05201   Mean   :0.8607   Mean   :0.8412  
##  3rd Qu.:0.5273   3rd Qu.:0.05754   3rd Qu.:0.8698   3rd Qu.:0.8637  
##  Max.   :0.6470   Max.   :0.13834   Max.   :0.9602   Max.   :0.9574  
##    cg05724065       cg11091262       cg09809672       cg11191337    
##  Min.   :0.4425   Min.   :0.5998   Min.   :0.2097   Min.   :0.5813  
##  1st Qu.:0.6296   1st Qu.:0.7612   1st Qu.:0.3798   1st Qu.:0.7602  
##  Median :0.6631   Median :0.7889   Median :0.4350   Median :0.7813  
##  Mean   :0.6579   Mean   :0.7869   Mean   :0.4327   Mean   :0.7790  
##  3rd Qu.:0.6916   3rd Qu.:0.8164   3rd Qu.:0.4832   3rd Qu.:0.8009  
##  Max.   :0.8539   Max.   :0.9646   Max.   :0.7617   Max.   :0.8941  
##    cg10092957        cg04183425       cg02228185       cg22580512    
##  Min.   :0.05526   Min.   :0.6014   Min.   :0.3342   Min.   :0.3187  
##  1st Qu.:0.10492   1st Qu.:0.7429   1st Qu.:0.6532   1st Qu.:0.6547  
##  Median :0.12311   Median :0.7732   Median :0.7279   Median :0.6976  
##  Mean   :0.12535   Mean   :0.7708   Mean   :0.7161   Mean   :0.6944  
##  3rd Qu.:0.13968   3rd Qu.:0.7991   3rd Qu.:0.7937   3rd Qu.:0.7421  
##  Max.   :0.29182   Max.   :0.8971   Max.   :0.9273   Max.   :0.8687  
##    cg06414605        cg16761581        cg07422345        cg07737781    
##  Min.   :0.03373   Min.   :0.02050   Min.   :0.03564   Min.   :0.1045  
##  1st Qu.:0.04758   1st Qu.:0.03815   1st Qu.:0.05589   1st Qu.:0.1683  
##  Median :0.05256   Median :0.04855   Median :0.06267   Median :0.1923  
##  Mean   :0.05319   Mean   :0.05022   Mean   :0.06456   Mean   :0.1955  
##  3rd Qu.:0.05783   3rd Qu.:0.06016   3rd Qu.:0.07032   3rd Qu.:0.2182  
##  Max.   :0.08362   Max.   :0.13078   Max.   :0.24595   Max.   :0.4361  
##    cg20645058       cg19904463        cg04474832       cg15804973    
##  Min.   :0.6959   Min.   :0.03365   Min.   :0.1556   Min.   :0.1418  
##  1st Qu.:0.8273   1st Qu.:0.07166   1st Qu.:0.2435   1st Qu.:0.2820  
##  Median :0.8413   Median :0.07854   Median :0.2648   Median :0.3337  
##  Mean   :0.8416   Mean   :0.07962   Mean   :0.2612   Mean   :0.3321  
##  3rd Qu.:0.8558   3rd Qu.:0.08565   3rd Qu.:0.2806   3rd Qu.:0.3833  
##  Max.   :0.9372   Max.   :0.18379   Max.   :0.3387   Max.   :0.5287  
##    cg01173291        cg22484793       cg17471102       cg04865575     
##  Min.   :0.01807   Min.   :0.6469   Min.   :0.2375   Min.   :0.01396  
##  1st Qu.:0.05613   1st Qu.:0.7155   1st Qu.:0.4747   1st Qu.:0.03348  
##  Median :0.06275   Median :0.7343   Median :0.5156   Median :0.03661  
##  Mean   :0.06333   Mean   :0.7336   Mean   :0.5104   Mean   :0.03765  
##  3rd Qu.:0.07045   3rd Qu.:0.7515   3rd Qu.:0.5481   3rd Qu.:0.04063  
##  Max.   :0.10261   Max.   :0.8474   Max.   :0.7356   Max.   :0.08213  
##    cg13931228        cg06493994       cg08090640       cg21801378     
##  Min.   :0.03480   Min.   :0.1221   Min.   :0.1797   Min.   :0.01947  
##  1st Qu.:0.09466   1st Qu.:0.2213   1st Qu.:0.4237   1st Qu.:0.06815  
##  Median :0.10839   Median :0.2410   Median :0.4559   Median :0.08571  
##  Mean   :0.11071   Mean   :0.2420   Mean   :0.4510   Mean   :0.09082  
##  3rd Qu.:0.12467   3rd Qu.:0.2605   3rd Qu.:0.4836   3rd Qu.:0.10849  
##  Max.   :0.19319   Max.   :0.3888   Max.   :0.6083   Max.   :0.29011  
##    cg02797569       cg03996822       cg18267374       cg15037004    
##  Min.   :0.6644   Min.   :0.2506   Min.   :0.0803   Min.   :0.3129  
##  1st Qu.:0.7553   1st Qu.:0.4115   1st Qu.:0.1309   1st Qu.:0.5019  
##  Median :0.7750   Median :0.4619   Median :0.1482   Median :0.5397  
##  Mean   :0.7732   Mean   :0.4596   Mean   :0.1508   Mean   :0.5397  
##  3rd Qu.:0.7932   3rd Qu.:0.5063   3rd Qu.:0.1666   3rd Qu.:0.5836  
##  Max.   :0.8856   Max.   :0.7201   Max.   :0.3073   Max.   :0.7122  
##    cg01511567        cg09817114     
##  Min.   :0.08287   Min.   :0.02313  
##  1st Qu.:0.16611   1st Qu.:0.05332  
##  Median :0.19103   Median :0.05807  
##  Mean   :0.19003   Mean   :0.05881  
##  3rd Qu.:0.21399   3rd Qu.:0.06360  
##  Max.   :0.29753   Max.   :0.09274

# limits of all CpG values
range(apply(mdna[, vec_probe], 2, min), 
      apply(mdna[, vec_probe], 2, max))

## [1] 0.002807284 0.993269100

# try: 

plot(apply(mdna[, vec_probe], 2, min), ylim = c(0, 1), pch = 16, col = 2, type = "l")
lines(apply(mdna[, vec_probe], 2, max), pch = 16, col = 3)

# the same thing but easier to look at
as.data.table(mdna) %>% 
  melt(measure.vars = vec_probe, variable.name = "id_probe", value.name = "methylation") %>% 
  .[, `:=` (m_min = min(methylation), m_max = max(methylation)), by = id_probe] %>% 
  .[, m_diff := m_max - m_min] %>%
  .[, .(id_probe, m_min, m_max, m_diff)] %>% 
  unique %>% 
  ggplot(aes(x = reorder(id_probe, m_diff))) +
  geom_segment(aes(y = m_min, yend = m_max, xend = reorder(id_probe, m_diff))) +
  geom_point(aes(y = m_min)) +
  geom_point(aes(y = m_max)) +
  labs(x = "", y = "Methylation") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

The values are given as methylated percentage.

Question: Why should we transform the CpG values?

# function for the log-odds transformation
logodds_trans <- function(p) log(p / (1 - p))

# transform the percentage values into log-odds
probe_logodds <- apply(mdna[, vec_probe], 2, logodds_trans)

# create data frame to work with
mdna_logodds <- data.frame(mdna[, 1:4], probe_logodds)

# alternatively, try with the methylation values as percentages
mdna_percent <- data.frame(mdna[, 1:4], mdna[, vec_probe]*100)

Henceforth we will be using the mdna_logodds data.

# CpG probes correlation with Age
vec_cor_values <- sapply(seq(vec_probe), function(x) cor(mdna_logodds$Age, mdna_logodds[, vec_probe[x]], method = "spearman"))

# look at the correlation values between each probe and Age
plot(vec_cor_values, ylim = c(-1, 1), ylab = "Spearman's Correlation")
abline(h = 0, lty = 2)

# density plot
plot(density(vec_cor_values), xlab = "Spearman's Correlation", xlim = c(-1, 1))

# save the CpG probes ordered by their absolute correlation value
vec_probe_cor <- vec_probe[order(abs(vec_cor_values), decreasing = TRUE)]

Question: What does having a high positive/negative correlation might indicate?

Questions

1. Main question: Can we predict one’s biological age?
2. Can we predict the age of an individual based on the methylation level of a given CpG probe?
3. Does gender pose an effect in this prediction?
4. Does ethnicity pose an effect in this prediction?

Exercise 1: Simple linear regression

For this exercise, we will study the relation between age and methylation levels from specific probes.

• Perform a simple linear regression to study the (linear) relationship between Age and a single specific marker (hint: use the function lm())
• For each case, what is the model’s formula (hint: coef() or coefficients())
• Model diagnosis
  • Look at the residuals (hint: resid() or residuals())
  • Assess their normality (hint: shapiro.test())
• Interpret the results
  • What are the models’ slopes?
  • What are the model predictions?
  • What are the predicted ages of someone with, for any inferred marker, methylation levels of 10%, 25%, 50%, 75%, and 90%? (hint: predict() with the newdata argument)

First probe

Linear model using the first CpG probe, cg17568996.

# CpG: cg17568996

# linear regression: formula as Y ~ x
fit1 <- lm(Age ~ cg17568996, data = mdna_logodds)

# look at the models' summary
summary(fit1)

## 
## Call:
## lm(formula = Age ~ cg17568996, data = mdna_logodds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.429  -9.449   0.610  10.725  38.191 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   49.434      3.618  13.665  < 2e-16 ***
## cg17568996    -7.277      1.781  -4.087 4.92e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.56 on 654 degrees of freedom
## Multiple R-squared:  0.0249, Adjusted R-squared:  0.02341 
## F-statistic:  16.7 on 1 and 654 DF,  p-value: 4.917e-05

Questions:

1. How do you interpret the Estimate values?
2. Which hypotheses tests are the p-values shown here referring to? What can we infer from the results?
3. What is the difference between the Multiple R-squared and Adjusted R-squared?

We can look at the model’s residuals distribution. Ideally, the distribution should be symmetric around zero (Normal distribution with \(\mu = 0\)). We can assess this visually, through the use of a density plot, histogram, violin plot, boxplot, or even a quantile-quantile plot (Q-Q plot).

# residual analysis
fit1$residuals

##            1            2            3            4            5            6 
##  -0.06529017  19.93009110   2.12637613  -1.64853927  -1.76916387  27.46043368 
##            7            8            9           10           11           12 
##   8.78665645  10.15608751   8.51566811  17.39883000  16.78307204 -19.90628754 
##           13           14           15           16           17           18 
##  15.61874132  -7.86584243  -4.89067693   6.62324274   2.67807288  14.27138696 
##           19           20           21           22           23           24 
##  24.80762830  -6.14360731  18.69478760   4.10029556  28.51311140  12.01075107 
##           25           26           27           28           29           30 
##  11.76707143   9.97469069  12.52975940  -7.13618154   6.23712849  10.87008250 
##           31           32           33           34           35           36 
##  14.40312830   1.50222691  19.61112712   6.63240748  10.66882976 -11.57921444 
##           37           38           39           40           41           42 
##  34.47579348  16.06932307  22.63642847  22.41715140  19.45518294   7.43948682 
##           43           44           45           46           47           48 
##  10.58439335   1.33637088  23.63891269  20.82560262  24.27034575   2.42402000 
##           49           50           51           52           53           54 
##  17.64037910  -7.38654054   7.55982620   2.32742907   7.71727358 -12.59470974 
##           55           56           57           58           59           60 
##  10.22839271  38.19064916  -7.87933934  11.33074480 -17.38465292  -9.69081642 
##           61           62           63           64           65           66 
##  10.58126428  21.02341370  -0.14603374  18.74279698   4.94372280   0.27705709 
##           67           68           69           70           71           72 
##  13.57278092   7.36224131   3.16397112  -6.61049816  11.31132867  13.86752390 
##           73           74           75           76           77           78 
##   6.13836013  21.45875724   9.09205314  18.62123441   2.27495136  24.28065228 
##           79           80           81           82           83           84 
##   9.11761965  -2.31382265  13.13067575   0.66873220  -1.47248128  22.50476048 
##           85           86           87           88           89           90 
##  20.60391447   4.49294505  11.80328427  10.71396418  16.79642590   6.22379503 
##           91           92           93           94           95           96 
## -11.31455281   7.00739651  -6.34805128  -0.55473152  20.18513095  11.90406339 
##           97           98           99          100          101          102 
##  14.04758966  -8.87332858   5.86421361  -0.89347518   4.47321518  15.91694548 
##          103          104          105          106          107          108 
##   6.95769459 -19.61288630   0.20711160   6.02152668  16.06605997  10.94066562 
##          109          110          111          112          113          114 
## -12.98931731 -12.55575549  -8.65979586  -6.29805793  -1.78635956 -17.55495862 
##          115          116          117          118          119          120 
##   4.85530060   9.80063428  -4.42390897  11.26224630 -10.04024592  -6.66696133 
##          121          122          123          124          125          126 
## -14.80114828 -20.74318966  -5.17856100  -0.69529957  16.71345490  -4.21232874 
##          127          128          129          130          131          132 
##   7.86892819  -2.25954677 -17.07983750  22.11935377   2.16205857 -13.38165924 
##          133          134          135          136          137          138 
##   7.27653738  13.41682097 -10.05056049  -9.34321008   0.42049319   5.88034771 
##          139          140          141          142          143          144 
## -16.34508130  -8.95518816   8.85852360 -16.02419610  16.85317677   5.78042911 
##          145          146          147          148          149          150 
##  -8.93797774 -14.36155528   3.78008353  -3.54415671 -16.74211030 -21.91964622 
##          151          152          153          154          155          156 
##  -9.13746187  -9.85267639  -8.40303835 -32.89336604  -2.07854270 -14.50501965 
##          157          158          159          160          161          162 
## -28.33929957 -18.49913413 -11.68758611  -5.18030930 -40.69450891 -11.59818967 
##          163          164          165          166          167          168 
##  -4.35317850  -8.91659300 -14.44801414 -30.84387177 -13.50672654 -17.76599194 
##          169          170          171          172          173          174 
## -16.50979300 -34.15985506 -27.86125963 -28.92639327 -16.85463102  13.75597504 
##          175          176          177          178          179          180 
##  23.68642048   5.13911973  13.99141715   5.57685270  21.27330682  18.16152036 
##          181          182          183          184          185          186 
##   7.39567525   3.38534774  -2.49432166   8.68473022  -6.40043814   0.95206297 
##          187          188          189          190          191          192 
##   7.86146466   9.17751281  11.33798594  11.40397737  10.75901932  14.57493317 
##          193          194          195          196          197          198 
##  -1.54236105  19.01815669  11.85332476  15.53000087  14.62020602  13.70228823 
##          199          200          201          202          203          204 
##   5.83634537  18.71781192  19.40329242  12.94603809  14.03197980  -1.97249418 
##          205          206          207          208          209          210 
##  -2.13182317   6.96768342  28.07522490   5.98424156   3.95682343  12.34514430 
##          211          212          213          214          215          216 
##  -1.05189809  15.84057046  16.08157060  19.86841152   2.26233610   4.85324701 
##          217          218          219          220          221          222 
##  -1.19408243  19.28053675   3.21188825   2.11037166  15.33041986   5.58401810 
##          223          224          225          226          227          228 
##  13.50986107   8.49532773   2.36843116   9.02762003  -5.66742232  10.32088312 
##          229          230          231          232          233          234 
##  20.35914471   5.50961122  -5.73216454   2.92071870   7.95463581   5.29879692 
##          235          236          237          238          239          240 
##  10.18726519   8.92236232  22.71339688  19.92866831  10.33826243   8.00991448 
##          241          242          243          244          245          246 
##  -2.98052071   2.24076538  18.26202778 -18.77229443 -20.25082580 -12.89616446 
##          247          248          249          250          251          252 
##  -1.99639393 -11.58319492 -16.97791093  -6.43026259 -25.02622073 -18.98550972 
##          253          254          255          256          257          258 
## -19.19771511 -11.81273531 -17.76483610   9.60645130 -11.56586495   1.13825374 
##          259          260          261          262          263          264 
## -17.61318042 -10.39457274   5.33553948 -16.07945191  -2.42709485  -6.64937690 
##          265          266          267          268          269          270 
## -20.66458863 -10.72423116 -14.91923661  -4.51815880 -10.88243249   0.20532649 
##          271          272          273          274          275          276 
## -23.03414269 -21.04323605  -2.50262035  -6.51966163 -18.41329890 -27.70394947 
##          277          278          279          280          281          282 
##  -9.66095270   7.67489310 -20.14229989  -2.21418813 -15.83569216 -23.47623534 
##          283          284          285          286          287          288 
##  -9.82454723   1.58604948 -15.16013448 -10.85754049  -0.69993418  -4.15333374 
##          289          290          291          292          293          294 
## -24.65649958  -7.44867093   0.80626769  -8.02482722  -5.65651168 -15.53178108 
##          295          296          297          298          299          300 
##   0.48770696  -0.51455489 -10.29319681 -18.66913175 -14.32528891  -9.45457654 
##          301          302          303          304          305          306 
##  -4.72668708  -0.21029935  -5.99314939 -13.60261220  -6.31703383 -29.93572191 
##          307          308          309          310          311          312 
## -17.13588314 -18.17191688  -0.90982139  -8.05663620  -0.39572376 -22.14016325 
##          313          314          315          316          317          318 
##   2.30133838 -16.70710552  -3.85730525  -1.98795760 -13.22351846 -17.01864401 
##          319          320          321          322          323          324 
##  -0.50136412 -16.32901452 -27.31620645   0.33263414 -10.26881170 -16.69554520 
##          325          326          327          328          329          330 
##  -4.01855548   0.98154021  -7.58859662 -26.95216231  -2.77983742  -1.91083366 
##          331          332          333          334          335          336 
##  -3.20581797  -8.73666152  -3.09628160 -12.76212764   0.14068693   2.94617346 
##          337          338          339          340          341          342 
##   0.77006797 -18.66162655  -1.07847605   4.18052434 -11.18532363  -5.09524921 
##          343          344          345          346          347          348 
## -10.32258839  -2.98417882  -3.99443361  -4.96201733 -10.29605859  -4.66770225 
##          349          350          351          352          353          354 
##  -3.68859963  -3.98338626  -0.95078527 -13.41992719 -32.70851243  -7.39373093 
##          355          356          357          358          359          360 
##  -9.44761304 -23.70554507  -3.96360551 -10.25872060  -8.96754078 -17.74196981 
##          361          362          363          364          365          366 
##   1.28456446  -2.79700895  -1.11573292   1.37265711  -0.70927397   2.72448907 
##          367          368          369          370          371          372 
##  -2.83626642  -3.54554523  -8.53681430  -6.27523760  -3.29748484 -11.16346441 
##          373          374          375          376          377          378 
##  -7.05885014   0.61810505 -17.78170376  -3.59856474  10.60188818 -15.78735047 
##          379          380          381          382          383          384 
## -15.37102591 -18.29128131  -2.91935785  -8.95271479 -10.69647773 -14.53188499 
##          385          386          387          388          389          390 
##  -9.10850029  -7.63379436   2.11419944   1.54044402   1.56535417  -8.23503550 
##          391          392          393          394          395          396 
##  -0.02637495   2.44701358   2.19173303  -3.60406090   2.45718478 -22.09564628 
##          397          398          399          400          401          402 
##  -6.10382218 -24.24372655   2.17720963  -5.33686740  -5.03847577 -21.40058566 
##          403          404          405          406          407          408 
## -11.25515228 -36.94273688  -1.60617318  -6.14502059   4.03366639  -0.10671465 
##          409          410          411          412          413          414 
##  -7.61541901  -8.91946807  -2.79991533  -6.28979493  -9.87404479   0.83352669 
##          415          416          417          418          419          420 
##  -5.60858091  -3.54461454   6.02282499   2.11301616  -1.69927216 -22.43811594 
##          421          422          423          424          425          426 
##  -7.34211425 -15.76725613 -13.72694913 -17.92819078  -9.01458937  15.14990632 
##          427          428          429          430          431          432 
##  16.04593596  12.66127013  13.32856623  -2.43108863   2.83112679  21.21336533 
##          433          434          435          436          437          438 
##  -8.79729979   4.87749181  21.57564645   4.86117802   3.20274339  22.27352371 
##          439          440          441          442          443          444 
##  13.67943212   2.04649685  15.00441466  -9.61961009  -3.25233703   3.56130148 
##          445          446          447          448          449          450 
##  -7.30947912  -6.14481494   4.83329643  -5.00001924  -2.72748193  27.59082476 
##          451          452          453          454          455          456 
##  13.47494202  24.64145899   5.98564849   4.92252587   0.73868019  26.21171283 
##          457          458          459          460          461          462 
##   2.16938463  17.51517172  -5.80010468  14.56202817  15.04707811 -16.31561237 
##          463          464          465          466          467          468 
##  27.29100510   0.60179948  -4.45731509   9.47797496  -4.67593066 -15.55670818 
##          469          470          471          472          473          474 
##   3.20530132   8.76660757  -0.46026879  18.33392903 -12.86216619  15.76320428 
##          475          476          477          478          479          480 
## -12.26245410   3.56887217 -12.19570978   7.78221980 -13.77630691  -2.87892513 
##          481          482          483          484          485          486 
##  -3.52670250   6.05330842  15.14451740  -6.60678138  15.82073469  16.49621399 
##          487          488          489          490          491          492 
##  -3.85595584  12.77868159  17.64551041  -7.92413987   3.11934323  15.35891716 
##          493          494          495          496          497          498 
##  -6.75443392   2.66574952  12.14558087  22.29805835   6.87779502  15.12468438 
##          499          500          501          502          503          504 
##  14.79084706   8.78740961 -16.47087313  -1.90361072 -10.74135952   2.88249866 
##          505          506          507          508          509          510 
## -25.21301631   6.13164902   4.08077895 -16.19795026   2.27170852  16.26357085 
##          511          512          513          514          515          516 
##  10.85699123 -16.78947876   4.07606404 -12.07728087  18.74566924  -7.74187473 
##          517          518          519          520          521          522 
##   2.88709166   8.87801321 -29.40744499   6.21356963   8.95700793  15.42405321 
##          523          524          525          526          527          528 
##   5.46888893  14.15969139   5.34671329 -13.32545410   1.42225663 -35.14131423 
##          529          530          531          532          533          534 
## -36.76558259 -16.62857656 -23.79335738 -17.18692446 -35.76409951 -15.53966331 
##          535          536          537          538          539          540 
##  -7.76652891 -26.67777876   5.02196443 -10.25441987 -18.24546129  -6.84888915 
##          541          542          543          544          545          546 
##   0.37486772   8.88144531  -6.51838558 -15.77461320  -6.28177327 -31.05955906 
##          547          548          549          550          551          552 
## -40.04774690 -21.02753613 -15.35973760 -19.40899468 -25.47533584 -35.88573313 
##          553          554          555          556          557          558 
##   5.59418383 -36.13682121 -11.54945892  -8.20102152 -45.42897833 -18.62932228 
##          559          560          561          562          563          564 
## -13.00138241 -28.85294263 -19.27870256  -6.66407664 -35.12160345 -17.15509921 
##          565          566          567          568          569          570 
## -29.41428494 -39.46512765 -24.03564050 -42.16317680 -23.51986965 -41.52408950 
##          571          572          573          574          575          576 
## -29.94048629 -11.26431797  16.51006503   3.93375743  22.50657215  26.46602603 
##          577          578          579          580          581          582 
##   0.80035142  -8.97963506  19.06264150   8.41752464  11.96801717  17.81420828 
##          583          584          585          586          587          588 
##  15.64018581  11.55943339  11.99416207  18.57641347  23.04793774  20.89330433 
##          589          590          591          592          593          594 
##  22.12269089  10.24750751  23.32259161  24.60485956  12.57581644  23.96717167 
##          595          596          597          598          599          600 
##  -4.46811252  10.38761672  13.35885367  20.86032013   4.26929811  18.81740243 
##          601          602          603          604          605          606 
##  25.07611013  21.29280039  23.51572439   0.56781198  22.54557900   8.09129673 
##          607          608          609          610          611          612 
##  20.81357716  18.06828251  21.32763913  31.03408508   5.46174001  21.22197340 
##          613          614          615          616          617          618 
##  20.57930599   8.61262100  13.50322456  -3.97168591  26.19282366  13.22513632 
##          619          620          621          622          623          624 
##   6.46302007  25.19560489  -7.94447077   2.32749787  13.44580088  31.67405027 
##          625          626          627          628          629          630 
##   3.53539417  12.97917755  21.25752696   1.74483900  -4.83145997   6.88508950 
##          631          632          633          634          635          636 
##  24.68181502  13.27969626   3.35116232  20.18868016  15.78796463   7.82366415 
##          637          638          639          640          641          642 
##  -0.26186194 -10.74637349   5.45946688   9.46781458  19.10028097   0.90377853 
##          643          644          645          646          647          648 
##   6.18647574  25.76861285  15.68191996   7.13914964   3.25805230  12.81251386 
##          649          650          651          652          653          654 
##  27.07850912   4.11909320  11.11214557  14.47344642   5.94834156   2.80386176 
##          655          656 
##  -2.23075251  11.03586152

# or
# resid(fit1)

# density plot
plot(density(fit1$residuals))

# Q-Q plot
qqnorm(fit1$residuals)
qqline(fit1$residuals)

# histogram
hist(fit1$residuals)

# diagnostic plots
par(mfrow = c(2,2))
plot(fit1)

par(mfrow = c(1,1))

Alternatively, we can apply formal hypothesis testing to infer on the residuals’ normality.

# proposal 1: Shapiro-Wilk normality test
shapiro.test(fit1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.99158, p-value = 0.0008728

# proposal 2: Lilliefors (Kolmogorov-Smirnov) normality test
nortest::lillie.test(fit1$residuals)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fit1$residuals
## D = 0.025678, p-value = 0.3674

Assuming this linear relation between age and the cg17568996 methylation values, what are the model fitted values (i.e., the prediction of the mean response value to the coefficient values used to formulate the model)?

pred_fit1 <- predict(fit1)


# look at the fitted results
boxplot(pred_fit1)

range(pred_fit1)

## [1] 56.35854 72.79467

plot(mdna_logodds$Age, pred_fit1, xlim = c(0, 100), ylim = c(0, 100))
abline(lm(pred_fit1 ~ mdna_logodds$Age))
abline(0, 1, lty = 2, col = "gray55")

What are the predicted ages of someone with cg17568996 methylation levels of 10%, 25%, 50%, 75%, and 90%?

# vector with percentages of interest
methylation_interest <- c(0.1, 0.25, 0.5, 0.75, 0.9)
# don't forget to transform the data!

# model predictions with 95% CI on estimates
predictions_fit1 <- predict(fit1, newdata = data.frame(cg17568996 = logodds_trans(methylation_interest)), interval = "confidence")

cbind(methylation = logodds_trans(methylation_interest), predictions_fit1) %>% 
  as.data.frame %>% 
  ggplot(aes(methylation, fit)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = "Log-odds Methylation", y = "fit1 predicted ages")

Top probe correlated with age

Perform the same analysis with the highest age-correlated probe, cg22736354.

# CpG probe: cg22736354

# linear regression
fit2 <- lm(Age ~ cg22736354, data = mdna_logodds)

# summary
summary(fit2)

## 
## Call:
## lm(formula = Age ~ cg22736354, data = mdna_logodds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.727  -7.964  -0.411   7.756  39.745 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  103.696      1.876   55.26   <2e-16 ***
## cg22736354    50.029      2.301   21.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.24 on 654 degrees of freedom
## Multiple R-squared:  0.4195, Adjusted R-squared:  0.4186 
## F-statistic: 472.6 on 1 and 654 DF,  p-value: < 2.2e-16

# residuals analysis
lillie.test(fit2$residuals)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fit2$residuals
## D = 0.031286, p-value = 0.1235

qqnorm(fit2$residuals)
qqline(fit2$residuals)

pred_fit2 <- fit2$fitted.values # this an alternative to the predict() of the fitted values

hist(pred_fit2)

# fitted results
plot(pred_fit2 ~ mdna_logodds$Age, xlim = c(0, 100), ylim = c(0, 100))
abline(lm(pred_fit2 ~ mdna_logodds$Age))
abline(0,1, lty = 2, col = "gray55")

# model predictions with 95% CI on estimates
predictions_fit2 <- predict(fit2, newdata = data.frame(cg22736354 = logodds_trans(methylation_interest)), interval = "confidence")

cbind(methylation = logodds_trans(methylation_interest), predictions_fit2) %>% 
  as.data.frame %>% 
  ggplot(aes(methylation, fit)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2) +
  scale_y_continuous(limits = c(-15, 300)) +
  labs(x = "Log-odds Methylation", y = "fit1 predicted ages")

# justification to the slope
plot(mdna_logodds$cg22736354, mdna_logodds$Age, ylim = c(0, 100))
abline(fit2)

Comparing results

From the results of the two linear models what do you see?

• Differences in the summary, R-squared, or coefficients?
• How about the residuals?

# try:

par(mfrow = c(1, 2))
plot(mdna_logodds$Age, predict(fit1), xlim = c(0, 100), ylim = c(0, 100), xlab = "Age, years", ylab = "Fitted age, fit1")
abline(lm(predict(fit1) ~ mdna_logodds$Age))
abline(0,1, lty = 2, col = "gray55")
plot(mdna_logodds$Age, predict(fit2), xlim = c(0, 100), ylim = c(0, 100), xlab = "Age, years", ylab = "Fitted age, fit2")
abline(lm(predict(fit2) ~ mdna_logodds$Age))
abline(0,1, lty = 2, col = "gray55")

par(mfrow = c(1, 1))

Exercise 2: Multiple linear regression

For this exercise, we will perform the same analysis for the cg22736354 probe, only this time the variables for gender and ethnicity will be included in the model.

• Perform a multiple linear regression
• Model diagnosis
• Compare against the previous model

2.1 Adding the coefficient for gender

# CpG: cg22736354


# define gender as a factor with "female" as reference
mdna_logodds$Gender <- relevel(factor(mdna_logodds$Gender), ref = "F")

# fit the model
fit2_g <- lm(Age ~ cg22736354 + Gender, data = mdna_logodds)

# summary
summary(fit2_g)

## 
## Call:
## lm(formula = Age ~ cg22736354 + Gender, data = mdna_logodds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.378  -7.872  -0.663   7.815  39.546 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 103.4908     1.8944  54.629   <2e-16 ***
## cg22736354   50.2005     2.3119  21.714   <2e-16 ***
## GenderM       0.7040     0.8819   0.798    0.425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.24 on 653 degrees of freedom
## Multiple R-squared:  0.4201, Adjusted R-squared:  0.4183 
## F-statistic: 236.5 on 2 and 653 DF,  p-value: < 2.2e-16

# assessing normality of residuals
lillie.test(fit2_g$residuals)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fit2_g$residuals
## D = 0.028325, p-value = 0.2274

qqnorm(fit2_g$residuals)
qqline(fit2_g$residuals)

# predictions for age in each gender
pred_g_female <- predict(fit2_g, newdata = data.frame(cg22736354 = mdna_logodds$cg22736354, Gender = "F"), type = "response", interval = "confidence")
pred_g_male <- predict(fit2_g, newdata = data.frame(cg22736354 = mdna_logodds$cg22736354, Gender = "M"), type = "response", interval = "confidence")

2.2 Adding the coefficient for ethnicity

# CpG: cg22736354

# define gender as a factor with the caucasian group as reference level
mdna_logodds$Ethnicity <- relevel(factor(mdna_logodds$Ethnicity), ref = "Caucasian")


# fit the model
fit2_e <- lm(Age ~ cg22736354 + Ethnicity, data = mdna_logodds)

# summary
summary(fit2_e)

## 
## Call:
## lm(formula = Age ~ cg22736354 + Ethnicity, data = mdna_logodds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.345  -6.996   0.222   6.762  39.535 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       101.6968     1.7343   58.64   <2e-16 ***
## cg22736354         43.2221     2.2039   19.61   <2e-16 ***
## EthnicityHispanic  -9.6886     0.8805  -11.00   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.33 on 653 degrees of freedom
## Multiple R-squared:  0.5103, Adjusted R-squared:  0.5088 
## F-statistic: 340.2 on 2 and 653 DF,  p-value: < 2.2e-16

# assessing normality of residuals
qqnorm(fit2_e$residuals)
qqline(fit2_e$residuals)

# predictions for age in europeans and mexicans
pred_e_cauc <- predict(fit2_e, newdata = data.frame(cg22736354 = mdna_logodds$cg22736354, Ethnicity = "Caucasian"), interval = "confidence")
pred_e_hisp <- predict(fit2_e, newdata = data.frame(cg22736354 = mdna_logodds$cg22736354, Ethnicity = "Hispanic"), interval = "confidence")

rbind(
  data.frame(Age = mdna_logodds$Age, pred_e_cauc, Ethnicity = "Caucasian"), 
  data.frame(Age = mdna_logodds$Age, pred_e_hisp, Ethnicity = "Hispanic")) %>% 
  ggplot(aes(Age, fit, col = Ethnicity)) +
  geom_point() +
  geom_smooth(method = "lm")

## `geom_smooth()` using formula 'y ~ x'

# try:

plot(mdna_logodds$Age, predict(fit2_e), xlim = c(0, 100), ylim = c(0, 100), col = mdna_logodds$Ethnicity, pch = 16)
abline(lm(predict(fit2_e) ~ mdna_logodds$Age))
abline(0,1, lty = 2)

plot(mdna_logodds$Age, predict(fit2_e), xlim = c(0, 100), ylim = c(0, 100), col = mdna_logodds$Ethnicity, pch = 16)
abline(0,1, lty = 2)
abline(lm(pred_e_cauc ~ mdna_logodds$Age), lwd = 3, col = 1)

## Warning in abline(lm(pred_e_cauc ~ mdna_logodds$Age), lwd = 3, col = 1): only
## using the first two of 6 regression coefficients

abline(lm(pred_e_hisp ~ mdna_logodds$Age), lwd = 3, col = 2)

## Warning in abline(lm(pred_e_hisp ~ mdna_logodds$Age), lwd = 3, col = 2): only
## using the first two of 6 regression coefficients

2.3 Interaction Interaction cg22736354:Ethnicity

• Try to add an interaction term between the two covariates used in the model
• Estimate the values for \(\beta_0\) and \(\beta_1\) on each population group

# adding ethnicity interaction
fit3 <- lm(Age ~ cg22736354:Ethnicity, data = mdna_logodds)

# summary
summary(fit3)

## 
## Call:
## lm(formula = Age ~ cg22736354:Ethnicity, data = mdna_logodds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.183  -7.300   0.228   6.683  45.839 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     99.474      1.803   55.16   <2e-16 ***
## cg22736354:EthnicityCaucasian   40.735      2.346   17.36   <2e-16 ***
## cg22736354:EthnicityHispanic    51.100      2.151   23.75   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.49 on 653 degrees of freedom
## Multiple R-squared:  0.4947, Adjusted R-squared:  0.4932 
## F-statistic: 319.7 on 2 and 653 DF,  p-value: < 2.2e-16

# predictions for age in the two groups
pred_fit3_cauc <- predict(fit3, newdata = data.frame(cg22736354 = mdna_logodds$cg22736354, Ethnicity = "Caucasian"))
pred_fit3_hisp <- predict(fit3, newdata = data.frame(cg22736354 = mdna_logodds$cg22736354, Ethnicity = "Hispanic"))

rbind(data.frame(fit = pred_fit3_cauc, Ethnicity = "Caucasian"), data.frame(fit = pred_fit3_hisp, Ethnicity = "Hispanic")) %>% 
  ggplot(aes(fit, fill = Ethnicity)) +
  geom_density(alpha = 0.5)

plot(mdna_logodds$Age, predict(fit3), xlim = c(0, 100), ylim = c(0, 100), col = mdna_logodds$Ethnicity, pch = 16)
abline(0,1, lty = 2)
abline(lm(pred_fit3_cauc ~ mdna_logodds$Age), lwd = 3, col = 1)
abline(lm(pred_fit3_hisp ~ mdna_logodds$Age), lwd = 3, col = 2)

Question: What is the difference between this model and two independent simple linear regression models for each group?

Extending the analysis

# run a simple linear model on every CpG probe available


# define empty data frame
simple_lm <- data.frame()


# run for loop
for(i in seq(vec_probe)) {

  # define formula to use
  fit_formula <- as.formula(paste0("Age ~ ", vec_probe[i]))
  # fit linear model
  fit <- lm(fit_formula, data = mdna_logodds)
  # model summary
  fit_summary <- summary(fit)

  # create a simple data frame with information about the model
  out <- data.frame(coef = vec_probe[i],
                    estimate = fit_summary$coefficients[2,1],
                    estimate_p = fit_summary$coefficients[2,4],
                    r_rquared = fit_summary$r.squared,
                    r_rquared_adj = fit_summary$adj.r.squared,
                    aic = AIC(fit),
                    bic = BIC(fit)
                    )

  # join rows of each model in output
  simple_lm <- rbind(simple_lm, out)
}

simple_lm

##          coef     estimate   estimate_p    r_rquared r_rquared_adj      aic
## 1  cg17568996  -7.27702442 4.916938e-05 2.490250e-02  0.0234115205 5379.829
## 2  cg02245378  11.15753778 5.936644e-05 2.436950e-02  0.0228777061 5380.188
## 3  cg14093936 -42.19132362 1.651705e-36 2.162624e-01  0.2150640696 5236.517
## 4  cg15642326 -15.30728848 7.618916e-08 4.326057e-02  0.0417976621 5367.361
## 5  cg24991452 -14.81362483 9.524716e-16 9.402093e-02  0.0926356380 5331.599
## 6  cg04581377  -4.37975783 4.445757e-02 6.161651e-03  0.0046420208 5392.318
## 7  cg01820374 -43.57284457 3.036497e-65 3.593396e-01  0.3583600162 5104.284
## 8  cg06898923   0.01823472 9.949354e-01 6.165529e-08 -0.0015289902 5396.372
## 9  cg11579069   7.52713317 3.397692e-02 6.856055e-03  0.0053374866 5391.859
## 10 cg19560758  43.43070533 1.775111e-51 2.943375e-01  0.2932584840 5167.679
## 11 cg17385448  -4.48862440 1.246486e-01 3.601632e-03  0.0020780868 5394.005
## 12 cg03271907   4.02597899 1.310724e-02 9.373653e-03  0.0078589335 5390.194
## 13 cg04081402 -26.11188624 1.602966e-29 1.769873e-01  0.1757288787 5268.594
## 14 cg26614073 -34.13126198 5.877815e-44 2.561205e-01  0.2549830285 5202.277
## 15 cg17895149   7.56682164 3.054310e-02 7.133513e-03  0.0056153687 5391.676
## 16 cg22736354  50.02905530 2.789930e-79 4.194977e-01  0.4186101052 5039.599
## 17 cg07549715   3.55529788 2.932023e-01 1.689250e-03  0.0001627808 5395.263
## 18 cg17592812  -2.45262988 4.331868e-01 9.395186e-04 -0.0005880968 5395.756
## 19 cg18463417   4.49154073 2.053417e-01 2.451417e-03  0.0009261132 5394.762
## 20 cg26832211   2.72617201 1.965226e-01 2.549200e-03  0.0010240459 5394.698
## 21 cg05724065 -34.42853644 7.959056e-44 2.554336e-01  0.2542951120 5202.883
## 22 cg11091262   3.00245825 1.769681e-01 2.785558e-03  0.0012607651 5394.542
## 23 cg09809672 -28.48343288 8.774189e-64 3.527336e-01  0.3517439412 5111.014
## 24 cg11191337  -7.28737222 1.193175e-02 9.625869e-03  0.0081115354 5390.027
## 25 cg10092957  -4.40347558 3.616496e-02 6.694106e-03  0.0051752901 5391.966
## 26 cg04183425   7.93264758 3.400991e-04 1.944955e-02  0.0179502346 5383.488
## 27 cg02228185 -16.59908511 1.508165e-57 3.237696e-01  0.3227356572 5139.731
## 28 cg22580512 -23.35950391 3.356067e-50 2.879890e-01  0.2869002826 5173.554
## 29 cg06414605  -1.82126409 6.255019e-01 3.644476e-04 -0.0011640471 5396.133
## 30 cg16761581  -2.44452399 1.461920e-01 3.225911e-03  0.0017017918 5394.253
## 31 cg07422345   9.24582115 3.581648e-04 1.930426e-02  0.0178047266 5383.585
## 32 cg07737781   7.27017131 1.782619e-03 1.482534e-02  0.0133189590 5386.574
## 33 cg20645058   0.24519783 9.377598e-01 9.330328e-06 -0.0015197074 5396.366
## 34 cg19904463   1.74476356 6.173792e-01 3.818219e-04 -0.0011466462 5396.122
## 35 cg04474832 -52.15043739 1.558379e-51 2.946174e-01  0.2935388471 5167.418
## 36 cg15804973 -26.76746200 9.793987e-68 3.704519e-01  0.3694892787 5092.806
## 37 cg01173291 -10.94372623 5.894338e-04 1.790816e-02  0.0164064856 5384.518
## 38 cg22484793   9.45610129 1.813329e-02 8.507175e-03  0.0069911313 5390.768
## 39 cg17471102 -37.02243020 1.508010e-57 3.237699e-01  0.3227358692 5139.731
## 40 cg04865575   2.80234747 3.644317e-01 1.257912e-03 -0.0002692161 5395.547
## 41 cg13931228 -30.35164193 1.840311e-37 2.214764e-01  0.2202860299 5232.139
## 42 cg06493994  50.29814226 3.148574e-61 3.410174e-01  0.3400097527 5122.782
## 43 cg08090640 -35.42772877 6.183052e-56 3.160713e-01  0.3150255842 5147.157
## 44 cg21801378  19.05031086 4.523630e-44 2.567133e-01  0.2555768180 5201.754
## 45 cg02797569  -7.82834728 1.975103e-02 8.280389e-03  0.0067639984 5390.918
## 46 cg03996822 -29.22401629 1.006238e-50 2.905973e-01  0.2895125622 5171.147
## 47 cg18267374  34.46292890 4.422488e-45 2.619573e-01  0.2608288174 5197.110
## 48 cg15037004 -29.89526944 1.857451e-44 2.587251e-01  0.2575916598 5199.977
## 49 cg01511567 -31.19460600 2.227967e-37 2.210236e-01  0.2198325070 5232.520
## 50 cg09817114   6.43151352 8.618342e-02 4.495227e-03  0.0029730487 5393.417
##         bic
## 1  5393.288
## 2  5393.646
## 3  5249.976
## 4  5380.820
## 5  5345.058
## 6  5405.776
## 7  5117.743
## 8  5409.831
## 9  5405.318
## 10 5181.137
## 11 5407.464
## 12 5403.653
## 13 5282.053
## 14 5215.736
## 15 5405.134
## 16 5053.058
## 17 5408.722
## 18 5409.214
## 19 5408.221
## 20 5408.156
## 21 5216.341
## 22 5408.001
## 23 5124.472
## 24 5403.486
## 25 5405.425
## 26 5396.946
## 27 5153.189
## 28 5187.013
## 29 5409.592
## 30 5407.711
## 31 5397.043
## 32 5400.032
## 33 5409.825
## 34 5409.580
## 35 5180.877
## 36 5106.265
## 37 5397.977
## 38 5404.226
## 39 5153.189
## 40 5409.005
## 41 5245.597
## 42 5136.241
## 43 5160.615
## 44 5215.213
## 45 5404.376
## 46 5184.605
## 47 5210.568
## 48 5213.435
## 49 5245.979
## 50 5406.875